Multi-frequency spatial frequency domain imaging: a depth-resolved optical scattering model to isolate scattering contrast in thin layers of skin

Abstract. Significance Current methods for wound healing assessment rely on visual inspection, which gives qualitative information. Optical methods allow for quantitative non-invasive measurements of optical properties relevant to wound healing. Aim Spatial frequency domain imaging (SFDI) measures the absorption and reduced scattering coefficients of tissue. Typically, SFDI assumes homogeneous tissue; however, layered structures are present in skin. We evaluate a multi-frequency approach to process SFDI data that estimates depth-specific scattering over differing penetration depths. Approach Multi-layer phantoms were manufactured to mimic wound healing scattering contrast in depth. An SFDI device imaged these phantoms and data were processed according to our multi-frequency approach. The depth sensitive data were then compared with a two-layer scattering model based on light fluence. Results The measured scattering from the phantoms changed with spatial frequency as our two-layer model predicted. The performance of two δ-P1 models solutions for SFDI was consistently better than the standard diffusion approximation. Conclusions We presented an approach to process SFDI data that returns depth-resolved scattering contrast. This method allows for the implementation of layered optical models that more accurately represent physiologic parameters in thin tissue structures as in wound healing.


Introduction
The skin is one of the largest organs of the human body, and it acts as an interface to the outside world and as a barrier to protect from external pathogens and other hazards.A wound is formed when skin is damaged or otherwise unable to perform its functions. 1,2The most common method used for evaluating the progression of the healing process is visual inspection by a trained physician.4][5] This method is fast, repeatable, and non-invasive, but the dependence on the expertise of the clinical personnel makes it a non-reliable way to obtain a consistent diagnosis. 6A more specific method to obtain information on the wound is biopsy, followed by staining and cell histology.This approach gives more structural information at the cellular level, but it is invasive, extremely localized and non-repeatable (especially on patients in-vivo), so it is only performed when there is already suspicion of malignant tissue.An ideal technique for wound assessment would combine the advantages of both visual inspection and histology, being non-invasive, repeatable in time, and specific, while presenting an objective measurement of parameters that can be employed as a metric for the evaluation of the efficiency of different regeneration treatments.
Optical techniques are non-invasive, fast methods for performing measurements on tissue using light.They satisfy most of the prerequisites that we want from an ideal wound healing assessment method: from the way that light interacts with tissue, it is possible to quantitatively measure concentrations of biological chromophores and indirectly obtain structural information on a microscopic scale, over a large field of view.Optical techniques are already widely utilized in clinical practice, where they are used for determining blood perfusion and blood oxygen saturation, [7][8][9] detecting cancer, 10,11 and optimizing photodynamic therapy, 12,13 among other applications. 14,15patial frequency domain imaging (SFDI) is an optical imaging technique that makes use of structured light illumination to quantitatively measure the absorption (μ a ) and reduced scattering (μ 0 s ) coefficients in tissue and separate their effects by making use of the frequency-specific response of the tissue. 16The absorption of light contains information about the concentration of molecules that are capable of absorbing photons (e.g., melanin, haemoglobin, etc.), whereas light scattering is dependent on the density and size distribution of the microscopic structures (e.g., collagen fibers, cell nuclei, etc.).These properties make SFDI a potentially valuable tool for the diagnosis of wound healing as it is able to measure objective parameters related to tissue function (e.g., blood concentration and oxygenation) and morphology (e.g., cellular proliferation and tissue remodeling) in a fast, non-invasive, and repeatable way.The technique, in its current state, has a few limitations, most notable of which is that the models of light transport, used to quantify the absorption and scattering parameters, presume a semi-infinite, optically homogeneous geometry.When looking at the practical application of wound assessment, it is evident that a homogeneous model is not an accurate representation of the underlaying physiology: biological tissue is quite heterogeneous and contains multiple thin layers with different properties.For the purpose of overcoming the limitations of homogeneous models we introduce a two-layer model composed of a thin layer, with a thickness in the range of hundreds of microns (which simulates the new tissue growth) and a thick layer, which can be approximated by a semi-infinite geometry and simulates the underlaying wound site.The introduction of the two-layer model alone, however, is not sufficient.To discriminate between the optical properties between the top and the bottom layer, we also need a dataset that contains information spanning multiple penetration depths of light into the tissue volume.It was shown in previous investigations that, when acquiring data over wavelengths of light spanning the visible and near infrared, each will have a different penetration depth depending on the optical properties of the tissue, with longer wavelengths (e.g., near-infrared) normally having longer penetration depths compared with shorter wavelengths (e.g., visible). 17owever, in the context of wound healing, this approach is not as applicable because the absorption contrast between layers is less pronounced.Scattering contrast, on the other hand, would contain distinct information about depth-specific structural changes in the wound site.SFDI also has a unique property of allowing for varying the penetration depth of the light patterns independently from the wavelength of light.This is made possible by means of the spatial frequency (f x ) of the illumination patterns. 18Sinusoidal patterns with higher frequencies have a lower penetration depth compared with patterns with lower frequencies.By acquiring multiple datasets with an increasing average spatial frequency, we are effectively sampling smaller and smaller volumes, which contain information in increasing proportion from the superficial layer, compared with the bottom layer.This multi-frequency approach is particularly suited to detecting variations in scattering contrast over thin superficial layer, as absorption contrast is discernible mostly at low f x (and higher penetration depth), whereas scattering has a larger influence at high f x (smaller penetration depth).
Belcastro, Jonasson, and Saager: Multi-frequency spatial frequency domain imaging. . .In this study, we acquire an SFDI dataset at several spatial frequency and then subdivide it in multiple overlapping sub-sets containing increasing spatial frequencies.These sub-sets, containing data with different penetration depth by means of different f x , constitute the basis of the two-layered model proposed in this study and are used to determine the layer specific scattering coefficients, which for the practical application of wound assessment can be correlated to new cell growth (re-epithelialization) and morphological changes within the underlying wound bed (remodeling).To test the validity of the model, SFDI data were acquired on silicone multi-layered optical phantoms with known properties and then processed according to the previously described multi-frequency approach.These tissue simulating multi-layer phantoms were designed to mimic the relative, layer-specific optical properties observed during reepithelialization of a wound.Using the individual properties of the two layers of the phantoms and their geometrical parameters, we compared the experimental data with the proposed twolayer model using three different analytical models of light transport.The root mean square percentage error (RMSPE) of the modeled data was then calculated as an evaluation metric to determine the performance of the models in different ranges of parameters.

Two-layer Scattering Model
Thin layered models of skin, based on the difference between the absorption in the epidermis and dermis layers, are already used in clinical applications to account for the different absorption levels of epidermis and dermis. 17,19Scattering-based multi-layer models, however, are less common.SFDI can be used to measure both absorption and scattering in a tissue, but for the purpose of this study, we consider differences in the scattering properties between the layers.As mentioned, in the context of wound healing, the difference in absorption between the layers is less pronounced, whereas the scattering coefficient can be subjected to large variations that give insight on the changes in tissue morphology due to the healing process.A two-layer scattering model based on the spectral difference in depth penetrance between the layers was already developed by our group to investigate morphologic changes within skin pigmentation and melanoma; however, that approach was optimized to isolate tissue layer thicknesses of ∼0.2 to 4 mm. 20n the present work, we aim to mimic the physiology of a healing wound, in which sub-cellular organelles of different sizes and new epithelial cells give different sources of scattering-based contrast.The end goal is to be able to measure this scattering contrast between thin tissue layers (0.05 to 0.5 mm).
The geometry of our two-layer model is represented in Fig. 1.The thin top layer has a finite thickness d and a scattering coefficient μ 0 S top .The bottom layer is a semi-infinite slab, with scattering coefficient μ 0 S bot .The index of refraction (n) is also assumed to be the same in both layers to avoid unwanted reflections and refractions due to index mismatch.There are three parameters in this scattering model: μ 0 S top , μ 0 S bot , and d.When performing normal SFDI measurements, we obtain a single measured scattering coefficient (μ 0 s meas ), which is assumed to derive from a single, homogeneous, and semi-infinite geometry.We want to define a modeled scattering coefficient (μ 0 s mod ) that models the behavior of μ 0 s meas and takes into consideration the relative contributions of the three parameters of the two-layer model.A first assumption that we make is that μ 0 s mod is a linear combination of the two free scattering parameters, μ 0 S top and μ 0 S bot , given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 4 ; 6 4 4 where αðdÞ is a weighting coefficient representing the layer-specific partial volume contribution of scattering, which is dependent on the thickness of the thin layer (d).

Fluence-Based Partial Volume Contribution (α)
Next, we need to find a metric to model αðdÞ as the relative contribution to light scattering from the top layer.Symmetrically, (1 − αðdÞ) represents the relative contribution to light scattering from the bottom layer.The metric that we chose to model αðdÞ is based on light fluence (ϕ).An advantage in using light fluence is that it can be calculated using different models that give an expression of fluence in three dimensions, so they contain the information about depth that we need to model light transport in multiple stacked layers.Because our technique uses planar illumination and wide-field imaging for detection, we assume the fluence to be constant relative to the field of view of a given pixel in the xy directions, simplifying the problem to one dimension in the z direction.
In our equation, we used the photon hitting density, 21 which for our source-detector configuration corresponds to the fluence squared (ϕ 2 ), because this is a co-located forward-adjoint model.Both our light source and detector (planar illumination and wide-field imaging, respectively) have the same probability distribution for a photon to reach a certain depth z (represented by ϕðzÞ).When multiplying the two functions for the source and detector, we end up with the same function squared.Because we are interested in the relative photon contribution from each layer, we define αðdÞ as the integral of the fluence in depth over the top thin layer, normalized to the integral of the total fluence, given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 4 ; 3 6 6

Influence of Spatial Frequency (f x )
To separate the optical properties of the two layers, we want to perform measurements that investigate different volumes of tissue.To do so, we make use of a specific property of SFDI: the penetration depth of the sinusoidal patterns is inversely proportional to their spatial frequency. 18,22This is possible because we are performing measurements of μ 0 S meas on a twolayered geometry, so changing the penetration depth of light will change the relative contribution of each layer to the overall scattering of light, as opposed to homogeneous models in which every measurement is supposed to return the same value, independently from the penetration depth of light.In Eq. ( 1), the layers with scattering coefficients μ 0 S top and μ 0 S bot are assumed to be homogeneous, so their values are independent from f x .The parameter α, instead, is based on light fluence, so it is also dependent on the spatial frequency, in addition to the thickness of the thin layer (d); it is given as

Analytical Models of Fluence (ϕ)
The last part needed to define our two-layer model is to choose a suitable model of light fluence (ϕ) that solves the radiative transport equation, specific to the spatial frequency domain.There are several options available to calculate a statistically reliable representation of ϕ: in general, Belcastro, Jonasson, and Saager: Multi-frequency spatial frequency domain imaging. . .from analytical models that make use of approximations to obtain simplified mathematical formulations of ϕ [23][24][25] to computer simulations using the Monte Carlo method that use realistic geometries and random propagation for millions of photons. 26,27In this study, we selected and evaluated the performance of two analytical models of light fluence.This decision was made because of the ease of implementation and much faster computation times compared with simulation-based models.In addition, normal Monte Carlo methods cannot currently perform direct simulations in the spatial frequency domain (SFD), so the spatial frequency dependency is obtained indirectly by means of the Hankel transform, 28 which is only applicable in certain conditions, at the cost of losing all spatial information.

Standard diffusion approximation
The first analytical model that we investigated is the standard diffusion approximation (SDA). 23he SDA is a model already used in many medical applications that works best in regimes in which scattering is predominant compared with absorption (μ 0 s ≫ μ 0 a ). 16For this reason, the SDA is more effective for light in the infrared spectrum because of the low absorption from biological tissue.The formula of light fluence (normalized to the incident power P 0 ) in the SFD for the diffuse approximation is given in the following equation: where A and C are constants derived by the choice of an appropriate boundary condition and R is a constant dependent on the effective reflection coefficient, which are given as ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 7 ; 4 5 9 For the full derivation and explanation of the coefficients, we refer to Cuccia et al. 16

δ-P1 approximation
Because our SFDI instrumentation uses a light source that includes the visible spectrum and we are trying to detect photons that have short pathlengths, the SDA becomes less accurate as we are approaching the limits at which the model is still reliable.For this reason, we also considered a second analytical model.The δ-P1 approximation is a diffusion-based model that introduces a correction factor to account for photons that are non-diffuse.This correction factor allows for modeling fluence more accurately in those conditions in which the SDA becomes less reliable (e.g., for short distance photon propagation). 24,29The original δ-P1 model did not include dependency on the spatial frequency f x , 29 which was only introduced later, in the context of a doctoral thesis. 30The formula of light fluence (normalized to the incident power P 0 ) in the spatial frequency domain for the δ-P1 model is given in Eq. ( 6).For a complete derivation and explanation of the coefficients, we refer to the doctoral thesis of Seo. 30 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 7 ; 2 2 8 where Þμ s as described in Carp et al. 29 and the constants C Ã and z b are obtained by the application of the boundary condition, as described in the appendix of Seo's doctoral thesis and are given as 30 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 1 1 7 ; 1 2 6

Modified δ-P1 approximation
In this study, we also derived our own spatial frequency-dependent δ-P1 equation based on the original δ-P1 equation 29 and applied the same modification and boundary conditions made by Cuccia et al. in their derivation of the SDA in the spatial frequency domain. 16We call this derivation modified δ-P1, or mod-δ-P1.The equation of light fluence for the mod-δ-P1 approximation is shown in the following equation: where the coefficients A 0 and C 0 are obtained by solving Boltzmann differential equation with the appropriate boundary conditions as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 4 ; 4 0 1 For a full derivation and explanation of the coefficients, we refer to Appendix A.
A comparison between the SDA and the two δ-P1 models of fluence are given in Fig. 2, which shows the values of ϕðzÞ, normalized to the total integral of ϕ over z.The models were calculated for optical properties μ a ¼ 0.05 mm −1 and μ 0 s ¼ 5 mm −1 .Figure 2 also shows how the spatial frequency changes the distribution of the photons, with higher values of f x having photons more concentrated near the surface.

Thin Phantoms Fabrication
To test the validity of the models, we manufactured silicone multi-layered optical phantoms with controlled scattering properties and thicknesses.These multi-layer phantoms were used to acquire experimental data to which the models were compared.India ink was used as an absorber because of its relatively flat absorption spectrum across the entire measured range.Because this study is focused on developing a model of light scattering, we want to minimize the spectral features generated by absorbers, which could alter the scattering spectrum in the case of cross-talk during the processing of the data.The ink was mixed in a sufficiently large volume of unpolymerized silicone with a concentration of 20 mg∕100 ml, before splitting it into different batches of phantoms, so that the absorption would remain consistent across all manufactured phantoms.Given the ink concentration used, we expected a flat absorption value of about 0.0167 mm −1 , 31,32 which was also independently measured using optical methods (3.2.1), to account for imprecisions in the fabrication process and the reliability of the used phantom recipe.The silicone was divided in two to fabricate two different batches of phantoms with different scattering particles.The scatterers were mixed in the silicone curing agent and sonicated to break down any aggregation of particles that would change their size distribution.Then the curing agent was mixed with the unpolymerized silicone in a ratio of 1:5, and a vacuum chamber was used to remove any air bubbles formed in the process. 31From each batch, five thin phantoms, with a thickness approximately in the range 0.1 to 1 mm, were obtained using a syringe to measure small volumes of liquid silicone (0.5 to 3 ml), poured onto a petri dish of 5 cm in diameter, and left to spread by resting on a level surface.The thickness of each phantom was then measured using a micrometer, taking multiple readings in different spots (N ¼ 10).A summary of the value of thickness (mean and standard deviation) is given in Table 1.The rest of the silicone was poured into a container to obtain a thick homogeneous phantom of about 20 to 30 mm in thickness.The scattering agent used in the first batch was aluminum oxide (Al 2 O 3 ) to emulate the underlying wound morphology, with an average particle size of 5 μm, in a concentration of 3 g∕100 ml.
The scattering coefficient measured at 650 nm is expected to be ∼0.8 mm −1 , with a scattering slope of ∼0.15. 33The scattering agent used in the second batch was titanium oxide (TiO 2 ) to emulate cellular proliferation, with an average particle size of 200 nm, in a concentration of 130 mg∕100 ml.The scattering coefficient measured at 650 nm is expected to be ∼1.78mm −1 , 31 with a scattering slope of ∼1.2. 332 Data Acquisition

Phantoms characterization
First, the optical properties of the two batches were measured on the 30 mm thick homogeneous phantoms using a handheld SFDI imager, capable of getting spatial information at five spectral bands in the visible range (450 to 630 nm). 34Because there is no spatial heterogeneity in the phantoms' optical properties, a rectangle was drawn on the center of the SFDI images, and the spatial average was calculated.The same measurements were also independently repeated for validation using spatial frequency domain spectroscopy (SFDS), which is a probe-based system capable of doing measurements with high spectral resolution in the range 400 to 1000 nm on a single point in space. 35Both SFDI and SFDS measurements were calibrated with respect to the same reference homogeneous phantom of known optical properties.The SFDS measurements were taken in three different spots for each phantom, but no significant changes were detected, which shows the homogeneous composition of the targets.The scattering coefficients at 650 nm are ∼2.07 mm −1 for TiO 2 and 0.8 mm −1 for Al 2 O 3 , and the absorption coefficients of both were within 0.016 to 0.026 mm −1 .Eventual differences from the values expected from the phantom recipes in Sec.3.1 are to be attributed to imprecisions in the fabrication process.The μ a and μ 0 s of  the two batches are reported in Table 2 at three spectral bands (458, 536, and 626 nm), with the respective ratio between the μ 0 s of the top and bottom layers (scattering contrast).

Multi-frequencies measurements
After having characterized the two batches, data were acquired on multi-layered phantoms using the SFDI imager previously described.The two-layered phantoms were obtained by placing each of the thin phantoms over the thick homogeneous phantom from the opposite batch to have two layers with distinct scattering properties.To improve adhesion between the layers and to reduce the refraction index mismatch due to the presence of air bubbles, ultrasound gel was spread between the phantoms and was squeezed out as much as possible.The Al 2 O 3 thick phantom was used as the bottom layer, and the TiO 2 thin phantoms were the top layer.A total of 11 spatial frequencies were acquired for each dataset, from 0 mm −1 (planar illumination) to 0.5 mm −1 , in incremental steps of 0.05 mm −1 .The two-layered phantom measurements were calibrated to the same reference used in Sec.3.2.1 on homogeneous phantoms, minimizing the introduction of calibration errors.

Data Processing
Our new multi-frequencies processing approach of subdividing the dataset into several smaller, partially overlapping sub-sets containing four frequencies each, with increasing values of f x , as seen in Fig. 3.Each sub-set is individually processed to obtain μ a and μ 0 s values, which have different penetration depths, as defined in Sec. 2. The average f x of each subset is used in the fluence calculation, denoted as hf x i in the text.The SFDI data were processed assuming a homogeneous volume using the procedure described in Cuccia et al.: 16 the amplitude of the modulated patterns was extracted and a homogeneous reference phantom with known optical properties was used for calibration to obtain the measurements of diffuse reflectance (R d ).
A white Monte Carlo model of R d is then used in an iterative algorithm that changes the input parameters ðμ a ; μ 0 s Þ according to an optimization strategy, until the R d of the model matches the R d measured on the target, giving as a result the ðμ a ; μ 0 s Þ pair of the measured target.s of these two layers are assumed to be constant across the spatial frequencies and were obtained by measuring it on the thick homogeneous phantoms using the same SFDI system.

Fluence-Based Models
Using the values of μ 0 s top and μ 0 s bot measured in Sec.3.2.1 and the phantom thickness (d) reported in Table 1, the expected values of μ 0 s ðf x Þ of the two-layered phantoms were modelled using Eqs.( 1) and ( 3).The μ a and μ 0 s measured on the layered phantoms, as described in Sec.3.2.2, were used to calculate the fluence ϕ used in Eq. ( 3) with the three models presented in Sec.2.3.The procedure is summarized in the flowchart in Fig. 5.In Fig. 6, a comparison between the modeled μ 0 s values and the experimental measurement on phantoms at 536 nm is shown for both SDA and the two δ-P1 models.

Models' Performance
To obtain an objective evaluation of the accuracy of the two-layered models of μ 0 s and compare their performance, the RMSPE was calculated with respect to the measurements described in Sec.4.1 for each of the models and each of the phantom thicknesses using the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 7 ; 1 5 7 where N is the number of average spatial frequencies contained in the dataset (N ¼ 8), μ 0 s meas is the measured scattering coefficient, and μ 0 s mod is the scattering coefficient modeled using Eq. ( 1).The chart in Fig. 7 summarizes the RMSPE of three datasets having different scattering contrasts, defined as the ratio μ 0 s top ∕μ 0 s bot .

Discussion
By observing the μ 0 s meas data in Fig. 4, we can make a few observations.First, the value of the μ 0 s coefficients measured on the two-layered phantoms are contained between μ 0 s top and μ 0 s bot , which is consistent with the hypothesis made in Sec. 2 of μ 0 s mod being a combination of the other two scattering coefficients.Second, with an increasing thickness of the thin layer, the μ 0 s meas of the two-layered phantoms becomes closer to μ 0 s top .This is to be expected because the relative contribution of the top layer to the overall scattering becomes greater the thicker the layer is.Similarly, with the increase of the average f x , the μ 0 s meas of the two-layered phantoms becomes closer to μ 0 s top .This behavior is also in line with the SFDI property shown in Sec.2.2, where higher spatial frequencies have less penetration depth, so the relative contribution from the top layer becomes higher.A last consideration that can be done is that the dependency of the twolayered μ 0 s over hf x i seems to be non-linear in nature and tends to "saturate" for high enough values of hf x i or layer thickness.This saturation of μ 0 s can be explained by assuming that the penetration depth of light (δ) becomes smaller than the thickness of the top layer, so we are effectively measuring a single layer instead of two.Looking at the two-layer scattering models in Fig. 6, it is possible to determine where the underlying fluence models are the most accurate and what their limitations are.The SDA model [Fig.6(a)] seems to be more accurate for larger thicknesses, but it tends to underestimate the contribution of the top layer, especially at low hf x i.This is expected, as SDA assumes that the photons have long enough pathlengths to be considered diffuse (distance from the light source ≫1∕μ t ).This behaviour is also evident from the fluence plots in Fig. 2 that show a larger concentration of photons deeper in the tissue compared to the other models.The original δ-P1 model [Fig.6(b)] seems to have an accuracy similar to the SDA in the extreme cases (d < 0.1 mm, d > 1 mm), but it is more accurate in the middle range (0.3 mm < d < 0.7 mm) and does not overestimate the contribution of the bottom layer at low f x , thanks to the additional correction term that models photons with a low number of scattering events/short pathlengths.The modified δ-P1 model [Fig.6(c)] is the most accurate for very thin layers (d < 0.2 mm), but it seems to over-estimate the contribution of the top layer and becomes less accurate for thickness d > 0.5 mm.In Fig. 7, we compare the performance of the models by looking at the RMSPE calculated for all thicknesses and three different spectral bands.Because the two scattering agents (TiO 2 and Al 2 O 3 ) have a large difference in the scattering slope, it is possible to choose bands far enough apart, so the ratio μ 0 s top ∕μ 0 s bot is different and we can do an evaluation with different contrasts between the top and bottom layer.In the three bands that we reported (458, 536, and 626 nm), the ratio is equal to approximately 4, 3.5, and 3, respectively.
A first consideration is that all three models have a similar performance for d > 1 mm, which might be an additional indication that we are exceeding the limits of the two-layers model and measuring only the top layer.Second, both δ-P1 models perform better than the SDA in most cases, with the original δ-P1 model performing better in the range 0.3 mm < d < 0.7 mm and our modified δ-P1 model performing better in the range d < 0.3 mm.Third, the RMSPE seems to improve in all models for a smaller μ 0 s top ∕μ 0 s bot ratio (i.e., less contrast between the layers), which is the opposite of what was initially expected.However, when considering that this two-layer approach is based on measurements that assume a homogeneous fluence distribution, it becomes unsurprising that errors in estimating the layer-specific scattering coefficients would increase the more we deviate from this assumption (i.e., for increasing difference in optical properties).
Nevertheless, these preliminary results remain extremely encouraging in the context of the target application (assessment of wound healing), in which we are interested in detecting a difference in μ 0 s in very thin layers of cells (0.1 to 0.2 mm).In particular, our modified δ-P1 model has an overall good performance (RMSPE < 10%), which becomes even better in the physiological range of interest (RMSPE < 5%).The objective of this initial investigation was to evaluate how closely homogenous fluence models could match SFDI measurement of layered media over differing penetration depths.
In Sec.2.3, we stated our motivations behind the selection of the analytical models of fluence analyzed in this study (i.e., computation speed and ease of implementation).6][27] This allows for future improvement by adopting a model that is capable of representing the fluence of light more accurately in different conditions and for different geometries and optical properties.
Multilayer fluence models in the spatial frequency domain do already exist 36 and can be utilized as a second stage optimization, using the results from the homogeneous models as initial estimates of layer thickness and layer specific optical properties.
The current work was aimed at modeling the depth-resolved multi-frequencies SFDI data in the most accurate way possible.However, the method still relies on pre-existent knowledge about the optical and geometrical parameters, which are the ones containing useful information for diagnostic purposes.For the method to be of practical use, an iterative inverse-solving algorithm will be implemented.The algorithm will allow for estimating the scattering parameters and layer thickness from the raw multi-frequencies SFDI measurements, providing the information most useful to aid clinicians in diagnosis.Furthermore, by combining the frequency-dependent depth estimation used in this work with the wavelength-dependent depth estimation from a previous work, 20 we will be able to offer a multi-purpose suite of optical tools for the analysis of various skin conditions, ranging from burn wounds to melanoma and everything else in-between.

Conclusions
We have presented an approach to processing SFDI data using different sub-sets of spatial frequencies to obtain datasets that have different depth information.We then made use of this depth-enhanced data to develop and validate a two-layer model of scattering for thin layers, aimed at mimicking the physiology of a healing wound.The model is based on the relative layer contribution to μ 0 s , calculated through an integral function of the light fluence.The performance of three analytical models of fluence was analyzed in the study, but the method itself is model agnostic and can be used with any model of fluence, leaving room for further improvement.The performance of the analytical model themselves looks promising, with the δ-P1 models overall working better than the SDA model.The proposed mod-δ-P1 model also has an excellent performance for very thin layers, which is especially interesting for the target application of this study, which is to measure layers of skin that are on the order of 0.1 to 0.2 mm.

Appendix A
The initial derivation follows the same procedure seen in the original δ-P1 model. 29By inserting the δ-P1 phase function in the Boltzmann transport equation, we obtain the following governing equations: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 1 1 4 ; 2 0 2 where . The source term that we introduce is of the kind E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 4 ; 1 5 9 qðr; zÞ where R s is the specular reflectance and q 0 ðrÞ is the irradiance of the light source.For planar illumination, we would have a constant value q 0 ðrÞ ¼ q 0 .This is where our approach differs from the original δ-P1 derivation, as we introduce a sinusoidally modulated light source, as seen in Cuccia et al., 16 which is given as We can then make the same considerations about the linearity of the medium, giving in response a sinusoidal fluence rate with no phase shift, given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 7 ; 7 1 2 ϕðzÞ ¼ ϕ 0 ðzÞ cosðk x x þ AÞ cosðk y y þ BÞ: (15) By introducing Eqs. ( 14) and ( 15) in Eq. ( 11), after the opportune simplifications, we obtain a 1-D differential equation in the depth dimension, given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 7 ; 6 6 3 where is the modified effective scattering coefficient, which is the term that introduces the spatial frequency dependence in the fluence equation, as seen in the standard diffuse approximation. 16rom here on, we can continue following the derivations steps found in the appendix of Carp et al., 29 solving the differential equation in a manner similar to the case of planar illumination on a semi-infinite geometry, but substituting the coefficient μ 2 eff with μ 02 eff .There are two boundary conditions required.First, because of the conservation of energy, the intensity of the diffuse light field must be zero in regions at a large distance from the source, given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 1 1 7 ; 5 3 0 The second boundary condition is given by the conservation of the diffuse flux component normal to the surface, which is given as where , where R 1 is the first moment of the Fresnel reflection coefficient in nonpolarized light.This implementation of R 0 is just an approximation of the correct term that makes use of the first two moments of the Fresnel coefficient.We chose to use this implementation because it gives a better approximation in proximity of the source, at the cost of reducing the model accuracy in the far field.Usually, the R 0 coefficient is obtained directly from a polynomial expression, as seen in Carp et al. 29 For the one-dimensional case, the equation of the second boundary condition is reduced to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 1 1 7 ; 3 4 9 The boundary conditions allow for solving Eq. ( 16) and finding the diffuse fluence rate, normalized to the incident power P 0 ð1 − R s Þ, with R s being the specular reflectance, which is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 7 ; 2 7 5 where E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 1 1 7 ; 2 2 6 Finally, by adding the collimated (non-diffuse) fluence contribution to Eq. ( 20), which is given as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 1 1 7 ; 1 4 9 ϕðz; which is in agreement with the solution presented in Eqs. ( 8) and ( 9).

Appendix B: Additional Figures
Additional data from the study is reported in this appendix and is shown in Figs.8-10.In the main text of the paper, the results from the two-layer phantom combinations were presented at 536 nm as a demonstration of the three different models' performance over increasing spatial frequency sets (Fig. 6).Here, Figs. 8 and 9 also show the same model performance relative to the measured bulk scattering results at 458 and 626 nm, which have a scattering contrast ratio between the layers of approximately 4 and 3, respectively.
A particular consideration could be given to the absorption coefficient measurements shown in Fig. 10.Given the manufacturing procedure of the phantoms, we expect a similar absorption in both the TiO 2 and Al 2 O 3 phantoms, which is reflected in the measurements at low hf x i.However, with increasing hf x i, we begin to see an increase in μ a in the multi-layered phantoms, especially the ones with the smallest thickness.We believe that this unexpected deviation is due to a combination of (1) the low sensitivity of absorption at high spatial frequencies (as previously shown in a study from Cuccia 16 ) and ( 2) the low signal to noise ratio (SNR) of the spatial frequency dependent reflectance at high spatial frequencies, as can be seen by the large variance of the measurements.This unsettling propagation of error in μ a ðhf x iÞ estimation may be a source of concern regarding the reliability of the proposed models and methods.We found that absorption contributes minimally, if not negligibly, to the fluence estimation at higher spatial frequencies and, subsequently, to this fluence-based approach to interpreting the partial volume contribution of layer-specific scattering.In this higher spatial frequency range, both hf x i and scattering become the dominant factors.Nevertheless, further investigation is warranted to determine the   Belcastro, Jonasson, and Saager: Multi-frequency spatial frequency domain imaging. . .exact source of this unexpected variation in μ a ðhf x iÞ as it could benefit future instrumentation design (i.e., impact of a better SNR, dynamic range, and calibration methods in measurements) and model development (i.e., more robust light transport models to reflect increasingly sub-diffusive absorption events).

Fig. 1
Fig. 1 2D drawing of the geometry of the two-layer scattering model (left).A thin layer of thickness d with scattering coefficient μ 0 S top is laid on a semi-infinite layer with scattering coefficient μ 0 S bot .Both layers are assumed to be homogeneous and extend infinitely in the xy plane.The light fluence ϕ (right) decreases exponentially from the surface in the z direction (depth).The weighting coefficient αðd Þ is defined as the partial contribution of fluence squared over the thin layer of thickness d.The contribution is obtained by integrating ϕ 2 up to d and normalizing to the total integral of ϕ 2 .

Fig. 2
Fig.2Comparison between fluence (ϕ) models using SDA (left), the original δ-P1 approximation (center), and our mod-δ-P1 derivation (right).The values of ϕ were normalized to their total integral over z.In each plot, the effect of the spatial frequency (f x ) on the fluence distribution is also shown.All models of ϕ in the figure were simulated on an homogeneous geometry with the same optical properties: μ a ¼ 0.05 mm −1 and μ 0 s ¼ 5 mm −1 .
the data measured in Sec.3.2.2 and processed assuming homogeneous tissue, as described in Sec.3.3, is shown in Fig.4.The figure contains data from a single spectral band (536 nm) because the measurements at each wavelength are independent from one another.The graph shows the scattering coefficients of the five two-layered phantoms and how they change with hf x i.The horizontal axis is the average of the spatial frequencies contained in each of the eight sub-sets.The μ 0 s of the top layer (TiO 2light orange line) and bottom layer (Al 2 O 3dark blue line) are included for comparison.The μ 0

Fig. 3
Fig.3Subdivision of the dataset containing 11 spatial frequencies (f x ) in 8 smaller, partially overlapped sub-sets, containing 4f x each.The sub-sets have different penetration depths (δ), which are estimated considering the average f x of the sub-sets (hf x i), with lower hf x i resulting in higher δ.

Fig. 5
Fig. 5 Summary of the procedure to acquire depth-sensitive data using multi-frequency SFDI and use it as input to the two-layer scattering model.SFDI images are acquired at several spatial frequencies and divided in sub-sets containing four frequencies each.From each sub-set, optical properties are estimated and used to calculate the homogeneous fluence rate ϕ.Values of ϕ are then used in combination with the model parameters to estimate the partial volume contributions α and obtain a model of μ 0 s ðhf x iÞ.

Fig. 4
Fig. 4 Scattering coefficients (μ 0 s ) of the five two-layered phantoms (colored markers), measured at 536 nm.The solid lines are the μ 0 s of the top layer (light orange) and bottom layer (dark blue).The data on the horizontal axis represent the average spatial frequency hf x i of the eight sub-sets into which the data were divided.

Fig. 6
Fig. 6 Comparison of the three fluence-based two-layer models (dashed lines) with the data (markers) measured at 536 nm.The top (light orange) and bottom (dark blue) solid lines are the scattering coefficients of the TiO 2 and Al 2 O 3 layers, respectively.The three models are (a) SDA, (b) the original δ-P1 approximation, and (c) our mod-δ-P1 approximation.

Fig. 7
Fig. 7 RMSPE calculated on the difference between the scattering coefficient of the three fluencebased two-layer models and the experimental measurements.Each figure has a different scattering contrast, defined as the ratio μ 0 s top ∕μ 0 s bot and shows the RMSPE of the three models for the five top layer thicknesses (in mm).The scattering ratios are approximately a ¼ 4, b ¼ 3.5, and c ¼ 3.

Table 1
Thickness of the thin TiO 2 phantoms, measured with a micrometer in 10 different spots to calculate the average and standard deviation.

Table 2
Optical properties of the two batches measured on 30 mm-thick phantoms, reported at three spectral bands.The first batch contains titanium oxide particles for scattering (row 1-2), and the second one contains aluminum oxide (row 3-4).Both batches contain India Ink for absorption.The ratio of the μ 0 s of the top layer to the μ 0 s of the bottom layer is reported on the bottom row.